arXivil 502.03926V 1 [math.NA] 13 Feb 2015 


Using cylindrical algebraic decomposition and local 
Fourier analysis to study numerical methods: two 

examples 

Stefan Takacs 
Faculty for Mathematics, 

Research Group Numerical Mathematics (Partial Differential Equations), 

TU Chemnitz, Germany 
Email: stefan.takacs@numa.uni-linz.ac.at 


Abstract —Local Fourier analysis is a strong and well- 
established tool for analyzing the convergence of numerical 
methods for partial differential equations. The key idea of local 
Fourier analysis is to represent the occurring functions in terms 
of a Fourier series and to use this representation to study 
certain properties of the particnlar numerical method, like the 
convergence rate or an error estimate. 

In the process of applying a local Fourier analysis, it is 
typically necessary to determine the supremum of a more or less 
complicated term with respect to all frequencies and, potentially, 
other variables. The problem of computing snch a snpremnm 
can be rewritten as a quantifier elimination problem, which can 
be solved with cylindrical algebraic decomposition, a well-known 
tool from symbolic computation. 

The combination of local Fonrier analysis and cylindrical 
algebraic decomposition is a machinery that can be applied to a 
wide class of problems. In the present paper, we will discuss two 
examples. The first example is to compute the convergence rate 
of a multigrid method. As second example we will see that the 
machinery can also be used to do something rather different: We 
will compare approximation error estimates for different kinds 
of discretizations. 

Index Terms —Mnltlgrld; Fonrier analysis; Cylindrical alge¬ 
braic decomposition 


I. Introduction 

In this paper, we want to give some examples where the 
combination of cylindrical algebraic decomposition (CAD), 
as a tool from symbolic computation, and local Eourier 
analysis (LEA) yield helpful results. LEA was introduced by 
A. Brandt, who proposed to use Eourier series to analyze 
multigrid methods, cf. m. Eor a detailed introduction into 
LEA, see, e.g., M- LEA provides a framework to determine 
sharp bounds for the convergence rates of multigrid methods 
and other iterative solvers for problems arising from partial 
differential equations. This is different to classical analysis, 
which typically yields qualitative statements only. So classical 
convergence proofs for multigrid solvers, cf. 0, show that 
the method is convergent and that the convergence rates are 
uniformly bounded away from 1 for all grid sizes, however 
there is no sharp, nor realistic bound for the convergence rate 
given. Besides the analysis of linear solvers, the idea of LEA 
can be carried over to other applications, like the computation 


of approximation error estimates or the computation of inverse 
inequalities. 

LEA can be justified rigorously only in special cases, e.g., on 
rectangular domains with uniform grids and periodic boundary 
conditions. However, results obtained with LEA can be carried 
over to more general cases, see, e.g. ||2l. In cases, where such a 
extension is not possible, it can be seen as heuristic approach. 

To compute the quantities of interest using LEA, typically 
one has to compute the supremum of a more or less com¬ 
plicated term. The key for involving symbolic algorithms is a 
proper reformulation of the problem of computing a supremum 
as a quantifier elimination problem, which can be solved using 
a CAD algorithm, cf. 0. Understanding the combination of 
LEA and CAD as a machinery for analyzing a numerical 
method, we apply this machinery in the present paper to two 
examples, keeping in mind that there are more. 

The first example is related to the classical idea of analyzing 
multigrid solvers. In Sec.|I^ we will introduce a classical finite 
element framework for the Laplace equation and analyze a 
standard Jacobi iteration for solving the discretized system. 
There, we will introduce the reader to the finite element 
method to keep the paper readable also for non-numerical 
analysts. In Sec. we will extend the analysis to be able 
to learn about convergence properties of a multigrid solver. 
The given example is rather simple (and could be solved also 
without use of CAD, just per hand). However, we refer to other 
examples, where the terms get much more complicated, which 
make symbolic tools more interesting, cf., e.g., 0 and fSl . 

is 


The second example, which will be discussed in Sec. IV 


a new result. It is given to show that the machinery of LEA can 
also be extended to analysis beyond analyzing the convergence 
of a multigrid solver. We will see that the method can also be 
used to develop approximation error estimates. Moreover, we 
will see that LEA can capture any kind of discretization. To 
keep it simple, we will stay in the one dimensional case, so 
the terms, that have to be resolved using CAD, are rather easy. 
We will provide supplementary material that covers also the 
extension to two dimensions. There, one can see that in this 
case the terms get much more complicated. 

This list of examples is not complete. So, CAD has already 



1.2 


been applied earlier in the analysis of (systems of) ordinary 
and partial differential-difference equations, El, where the 
necessary conditions for stability, asymptotic stability and 
well-posedness of the given systems were transformed into 
statements on polynomial inequalities using Fourier or Laplace 
transforms. 

11. Finite element method and a simple iteration 

SCHEME 

We start our analysis with a simple example, the Laplace 
equation. For a given function /, we are interested in finding 
a function u such that 

-u''{x) = f{x) ( 1 ) 

is satisfied for all x G 17 := (0,1) and, moreover, the boundary 
condition u(0) = u(l) = 0 holds. 

The standard way of solving this, is to introduce a varia¬ 
tional formulation. Let be the standard Sobolev space 

of weakly differentiable functions and iLo(17) C be 

the space of functions that moreover satisfy the boundary 
condition u(0) = u(l) = 0. Then, the strong formulation ([T]) 
can be rewritten in weak formulation as follows: Find u € 
V := Lfo(17) such that 

j u {x)v'{x)dx = j f{x)v{x)dx (2) 

Jq Jo. 

for all u G y, cf. standard literature on finite elements, like lO. 

For any finite dimensional subset Vk C V, we can introduce 
a discretized problem: Find Uk G 14 such that 

/ u'k{x)v'k{x)dx = / f{x)vk{x)dx (3) 

Jn Jn 

for all Vk G Vk- The approach to use the same space, 14, 
for both, Uk and Vk, is called the Galerkin principle. This 
guarantees that Uk is the orthogonal projection of the exact 
solution u G V into 14- 

The easiest way to set up the space 14 is to choose the 
Courant element: Here the domain H is subdivided into inter¬ 
vals (in one dimension) or into triangles (in two dimensions). 
We call these intervals or triangles elements. The space 14 
consists of all globally continuous functions that are linear on 
each element. 

Each function in 14 can be characterized just by prescribing 
its values on the end points of the intervals or at the vertices 
of the triangles, respectively - we call these points nodes. This 
fact can be used to construct a basis: The nodal basis of 14 
is the collection of all functions ^ G 14 that take the value 
1 on exactly one of the nodes and the value 0 on all of the 
other nodes. One such basis function is visualized in Fig. [T] 

Having this basis, we can represent the functions Uk and Vk 
in terms of the basis: 

N N 

Uk{x) = ^ Uk,i(pk,iix), Vk{x) = ^ Uk,i(pk,ii^)^ 

i=l i=l 

where the functions Uk and Vk can be represented by the 
coefficient vectors Uk := {uk,i)iLi Vk := {vk,i)fLi- 



Fig. 1. Basis functions of standard Courant element 

The variational equality ([^ can be rewritten in matrix-vector 
notation as follows: 

vlKkUk = vll^, ( 4 ) 

for all Vk G where Kk := 

and := (f^ fk(x)ipk,z(x)dx)fij^. As 0 is supposed to be 
satisfied for all V/., it can be rewritten as follows: Find Uj. such 
that 

KkUk=lk- (5) 

To obtain a good approximation, it is often necessary to refine 
the intervals (or triangles) used for the discretization of the 
partial differential equation. In this case both, the number of 
unknowns and the condition number of the matrix Kk, grow. 
However, Kk has a nice property: it symmetric and positive 
definite. 

A simple linear iteration scheme to solve a matrix-vector 
problem 0 for Kk being symmetric and positive definite, 
is the (damped) Jacobi iteration. Assuming to be some 
starting value, the iteration procedure is given by 

:= + T(diagA:fc)-^(/^ - Kki^^'’), 

where r > 0 is a given damping parameter. For r = 1, we 
obtain the standard Jacobi iteration. 

As a next step, we are interesting in analyzing the con¬ 
vergence of the Jacobi iteration scheme. So, using the exact 
solution := we obtain 

^pn+i) _ ^ (J _ r{dmgKk)~^Kk){u^^'’ - ul) 

and further 

\\uir^^^-y*k\\K < \\I-x{dmgKk)~^ KkllKju^r^ -u*k\\K^, 

where Sk ■= I — T{diagKk)~^Kk is called the iteration matrix 
and II • \\k^ is the vector norm ||ufe||if, := or 

the associated matrix norm. We have 

ll^felk. = \\Kl/^iI - r{dmgKk)-^Kk)K;^^\ 

1/2 

where || • || is the standard Euclidean norm. As K^J (/ — 
T[dmgKk)~^Kk)K^ is symmetric, obtain further 

ll^fclk, = - T{di^gKk)-^Kk)K-^'^) = p{Sk), 

where p(-) is the spectral radius. 

To determine the spectral radius, we use LEA: We compute 
the spectral radius of Sk explicitly for a special case. We 
assume to have 






• an infinitely large domain fl (this neglects all influence 
coming from the boundary of the domain), 

which is 

• discretized using an uniform (equidistant) grid. 

For simplicity, here, we restrict ourselves to the one dimen¬ 
sional case. However, LFA can also be worked out for two or 
more dimensions, cf. Go). 

For such an equidistant grid, we can compute the stiffness 
matrix explicitly: 

\ 

2 -1 
-1 2 -1 
-12-1 
-1 2 

J 

where is the grid size (length of the intervals). 

As next step, we define for any frequency 9 G [0, 27r)‘^ a 
vector of complex exponentials 

:= (</>k,j(9))jez := 

and observe that 

( 0 ) = (-e-"‘ + 2 - e®‘) 4 (9) ( 6 ) 

"-" 

Kk{ey.= 

is satisfied, i.e., that is an eigenvector of K/,. In the 

LFA world, the eigenvalue Kk{9) is also called the symbol of 
Kk. 

Based on the symbol of Kk, we can determine the sym¬ 
bol (eigenvalue) of the iteration matrix Sk- First note that 
diagiFfc = and therefore diagKk{9) = So, we obtain 

^(0) = l-ryFfc(0) 

= 1 -^(-e-®‘-h 2 -e®‘) = 1 -t( 1 -cos 6 I). ( 7 ) 

As we have mentioned above, we are interested in p{Sk)- This 
spectral radius can be expressed using the symbol: 

q{T) ■.= p{Sk) = sup |S'fc(6')| = sup |1 - r(l - cos6i)|. 

SG[0,27I-) eG[0.27i-) 

By substituting the variable 0 by c := cos 9, we can completely 
eliminate the occurrence of trigonometric functions and obtain 

q{T):= sup | 1 -t( 1 -c)|. 

-1<C<1 

By definition, the supremum is smallest upper bound, i.e., the 
smallest A such that 

V-i<c<i - A < 1 - r(l - c) < A. (8) 

To determine the smallest A satisfying ([^, we have to elim¬ 
inate the quantifiers, i.e. to solve a quantifier elimination 
problem. 

A quantifier elimination problem is the problem to And 
a quantifier free formula that is equivalent to a quantified 
formula: 




Fig. 2. Reduction of the high frequency modes as function of r 


Quantified formula: 

{f^l]x-i • ■ ■ {Qn)x„ ■ 5 ‘t'n; yi • ■ • ; 2/771)5 

where Qi G {3,V} and A is a finite boolean 
combination of polynomial inequalities 




Quantifier free formula: B{yi ... ,ym), 
where S is a finite boolean combination of poly¬ 
nomial inequalities. 

The solution of such a problem is possible using CAD, cf. a, 
a. By applying a CAD algorithm to we obtain 

(r < 0 A A > 1 - 2t) V (0 < T < 1 A A > 1) 

V (t > 1 A A >-l-h2r) (9) 

Here, the smallest A satisfying a is piecewise given by the 
terms 1 — 2t, 1 and -1-1- 2t. So, we obtain 

( 1 — 2 t for T < 0 
g(r) = \1 forO<T<l 

— 1 -|- 2t for 1 < T. 

We observe that there is no choice of r such that g(r) < 1. 
This reflects knowledge on the Jacobi iteration (which is also 
true for other simple linear iteration schemes): the convergence 
is not robust in the grid size hk, so the convergence rate 
cannot be bounded away from 1. (Although, we did not have 
an explicite dependence on the grid size hk, the fact that 
we have considered an unbounded domain fl is equivalent 
to considering an infinitely small gird size.) 

It is known by intuition that simple linear iteration schemes 
reduce high frequency error modes. This statement can be 
formally expressed using LFA: Here, we only consider 9 G 
[0,7r/2) U [37r/2,7r) or, equivalently, 0 < c < 1. In this case, 
we obtain using the same arguments as above 

9sm(t) = sup |l-r(l-c)| 

0<c<l 

Again, we can compute using CAD (or still per hand) that 

r 1 — 2r for r < 0 
ysM(r) = < 1-T for0<r<§ 

[ — 1 -I- 2t for § < T. 







-4>k (<2/) 


This function is visualized in Fig. We see that qgM takes 
its minimal value ^ for r = |. 


III. Analysis of a multigrid solver 


In the last section, we have seen that the Jacobi iteration re¬ 
duces the high frequency error modes. The idea of a multigrid 
method is to use the fact that low frequency error modes can 
be resolved well also on a coarse grid. So, we combine the 
Jacobi iteration (or any other simple linear iteration scheme) 
with a coarse grid correction, which reduces the low frequency 
error modes. 

We assume to have for fc = 1,2, 3,... a hierarchy of grid 
levels, where a grid level k is obtained from grid level k — 
1 by uniform refinement, i.e., in the case of one dimension: 
by subdividing each interval into two equally sized intervals. 
Starting from an iterate the next iterate of the 

multigrid method on grid level k is given by the following 
three steps: 

• Pre-Smoothing: Compute 

+ T(diag Kk)~^ - Kk • 

• Coarse-grid correction: 

- Compute the defect and restrict it to 

grid level k — 1: 

Tt-l •— ^k-l yPk —k J ■ 

- Solve the following coarse-grid problem approxima- 
tively: 

Kk-i^^\=rt\. (10) 


- Prolongate to the grid level k and add the result 
to the previous iterate: 


(m,2) 


(m.l) 


+ Pk-iP 


(m) 

k-1' 


• Post-Smoothing: Compute 

:= ^ 

As we have nested spaces, i.e., 14_i C I 4 , there is canonical 
embedding from 14_i into 14, which is chosen as prolonga¬ 
tion operator Pk-i- 

If the problem ( [TOl l is solved exactly, we obtain the two- 
grid method. In practice, the problem ( fTO] ) is approximatively 
solved by applying one step (V-cycle) or two steps (W-cycle) 
of the multigrid method, recursively. Only on the coarsest grid 
level, is solved exactly. 

For computing the convergence rate of the multigrid solver, 
we set up again the iteration matrix Qk, which is the product 
of the iteration matrix Sk of the damped Jacobi iteration, of 
the iteration matrix Ck of the coarse-grid correction and, once 
more, of the iteration matrix Sk of the damped Jacobi iteration: 


Gk — SkCkSk-, 



Fig. 3. Canonical embedding of 14—1 into Vj, 


and, as in the last section, 

Sk= I - T(diag Kk)~^Kk. 

As in the last section, we are interested in computing 

q{T) = \\Gk\\Kk = p{Gk)- 

To be able to determine the symbol of the iteration matrix 
Gk,'we have to take a closer look onto the prolongation opera¬ 
tor Pk-i first. We recall that there is an isomorphism between 
K^, the space of coefficient vectors, and the function space 
14 . So, for each coefficient vector = {(I>k,ji0))j^z, there 
is a function (pk{d, ■) G 14 ^ which is assigned to it: 

= '^(j)kj{0)pk,j{x). 

jGZ 

By definition, Pk-i is the canonical embedding operator, 
which is visualized in Fig. 

The next step is to represent the function <pk-i{‘^0,x) as a 
linear combination of functions on the fine grid. We observe, 
that this can be done using the ansatz 

(j)k-ii20,x) = A(j)k{9,x) -P B(j)k{0 -P 7r,x). 

It is sufficient to consider the nodes xj = jhk only. First 
we consider the even nodes X 2 j, which are also nodes of the 
coarse grid: 

(j>k-li2e,X2j) = A4>k{9,X2j) + B(j)k{9 -PTT,X2j). ( 11 ) 
As the {(pk,i)i€:i,^ form a nodal basis, o is equivalent to 
(l)k-i,j{29) = A4>k,2ji.(^) + B4>k,2j{0 PT^), 
eJ2»i = -f 


and, finally, 

1 = A + B . 

Now, we consider the odd nodes X 2 j+i, which do not occur 
on the coarse grid: 

4 > k - l { 20 , X2j + l ) = A4 > k { 0 , X2j + l ) PB ( l ) k { 0 + TT , X2j + l ). (12) 
As the {pk,i)i£^i„ form a nodal basis, is equivalent to 

2 i4’k-i,j{20) -P 0fe_ij_|_i(20)) 

= A(l)k^2j + l{0) + B(j)k,2j + l{0 + t ) 


where 


Ck = I- Pk-iK-\P^_^Kk 








Fig. 4. Coai-se-grid function (26», x) in black and the two components Fig. 5. Convergence rate of the multigrid solver as a function of r 

I (1 + cos{9))tf>k {9, x) and | (1 — cos{9))(f>i^{9 + tt, x) in gray. 


and 

i g2(j+i)eij ^ ^g(2i+i)ei ^g(2j+i)(e+7r)i 

and, finally, 

- (e-®‘ + e®‘) =A-B. 

cos{0) — 

We obtain A = ^(1 + cos(0)) and B = ^(1 — cos(0)), which 
can be observed also in Fig. This allows to introduce the 
symbol of the prolongation operator: 


^i{0) = I 


1 + cos(0) \ 

1 — cos(0) J ' 


Here, the symbol cannot be understood as eigenvalue any¬ 
more. However, for all 0 = [0, 27r), the prolongation operator 
Pk-i maps the linear span, spanned by 

(13) 

to the linear span, spanned by 


0^(6»-b7r), (14) 

and the restriction operator P^_i maps the linear span, 
spanned by ( |T4l i, to the linear span, spanned by ( [T3] l. 

Having this, we can set up the symbol for the two-grid 
operator Qk- We make use of the fact that the multiplication 
of Qk with a vector in the linear span, given by the basis ( [T4l i, 
maps into the same linear span. So, we have to set up the 
symbol of Gk with respect to the two dimensional basis d- 
The symbol of Sk has been a scalar in the last section. This 
means that every frequency was preserved by the action of Sk- 
If we represent the symbol of Sk with respect to the basis ( [T^ , 
we just obtain a diagonal symbol: 


Sk{0) 


Kio) 

KiO + n) 


where Skid) is as defined in Q. Exactly the same way, we 
obtain the symbol /Cfc(0) based on Kk{0), given in (|^. Using 
this, we can determine the symbol of Ck, 

Ckio)=1- 


where A* is the conjugate complex of . Consequently, the 
symbol of Gk is 

^kio) = Skio)Ckie)Skio). 


Here, the computation of Gki^) and of piGkid)) is straight¬ 
forward. We obtain: 

piGkid)) = |(r- l)^-br(3r-2)cos^(6»)|. 

As in the last section, we are again interested in computing 
the supremum 

qir) = piGk) = sup |(r - 1)^-b r(3T - 2) cos^(6»)|, 

ee [0.271-) 

where we again substitute cos 0 by c and obtain 

qir) = sup |(r — 1)^-j-T(3r — 2)c^|. 

cG[-l.l] 

Also here, we can resolve the supremum using a CAD 
algorithm (or, still, per hand) and obtain 

( 1 — 4t -I- 4t^ for r < 0 

qir) = < 1 — 2r -I- for 0 < r < | 

1 — 4t -I- 4t^ for I < T. 

This function is shown in Fig. We see that q takes its 

minimal value ^ for r = |. 

So far, all computations had been so easy such that it 
would have been possible to do them per hand. However, 
the methodology presented in this section can be carried over 
to more complex (and more interesting) problems. The first 
extension would be to consider two or more dimensions. Here, 
one could represent everything use a tensor-product structure, 
cf. m. Consequently, one has to deal with tuples of d 
frequencies for d dimensional spaces. Also in this case, the 
9i can be substituted by Ci := cos(0i) and solved as discussed 
in this session. However, the complexity of the expressions 
(particularly in terms of the polynomial degree) grows very 
fast if d is increased. 

Besides that, the presented methodology can be extended 
to non-standard problems. This is of practical use because the 
convergence analysis has to be worked out for each problem 
class, separately. Here, LEA can be of great help. 

One example where the presented approach has been ap¬ 
plied in this fashion was in a in a joint work with V. PillweiiiM 
cf. Q, H, where LEA and CAD have been used to compute 
convergence rates of a multigrid solver for a system of PDEs 
which characterizes the solution of an optimal control problem. 
There, not only the robustness of the convergence rates in 
the grid size hk, but also the robustness of the convergence 
rates in a regularization parameter, which is part of the 

'Research Institute for Symbolic Computation, Johannes Kepler University 
Linz, Austria 





problem description, was of interest and could be studied. 
The supplementary material, that came with the cited paper, 
is available in the wet^ The author wants to refer the reader, 
which is interested in analyzing multigrid convergence, to that 
material. 

In the following of the present paper, the author wants to 
draw the reader’s attention to another application of LFA that 
is also of interest in numerical analysis: the estimation of 
approximation error estimates. 

IV. Estimate the approximation error 

In this section, we are interested in comparing estimates of 
the approximation error 

inf \\u-Uk\\L^(Q) 

«fc6t4 

for different kinds of discretizations. One of the discretizations 
will be the Courant element, two more will be introduced 
below. Here and in what follows || • is the standard 

L^-norm, i.e., |l/||i 2 (n) :=J^Pix)dx 

One important approximation error estimate reads as fol¬ 
lows: 

inf ||u —< C'A^fc|u|^i(n) 

for all u G where Ca > 0 is a constant, hk is the grid 

size and |rt|rri(o) := ||w'||i 2 (Q). For classical discretizations, 
it is well-known that such an estimate exists. However, often 
there is no realistic bound for the constant Ca- So, it might 
be of interest to compute an realistic (not necessarily sharp) 
upper bound for the constant Ca for discretizations of interest. 

The approximation error can be bounded from above using 
an interpolation error ||m — where Hj. : -G 

Vk is an arbitrarily projection operator. So, it suffices to 
estimate 

||u — nfeu|||2(Q) < (15) 

for any projection operator H^. Using the following lemma, 
we show ([BJ for Hfc being the -orthogonal projection. 

Lemma 1: Let for all grid levels fc € N, the operator 
Hfe be the -orthogonal projection form into Vk- 

Assume that for all k the following quantitative estimate on 
two consecutive grids is satisfied: 

||(/— nfc)Mfc_i_i|||2(Q) < (16) 

for all Uk+i G Vk+i- Moreover, we assume to know qualita¬ 
tively that 

||(/- nfc)M||i 2 (o) ^ 0 for fc oo (17) 

for all u G Then the following estimate is satisfied: 

\\{i — nfc)M|||,2(Q) < 

for all u G 

Proof: The proof is based on a simple telescoping argu¬ 
ment. Due to ( [T7| ), for any e > 0 there is some AT > 0 such that 
11(7 — n^jwllj;^ 2 (f 2 ) < e|u|//i(f 2 ). Now, we obtain due to the 


triangular inequality, ([T^ and the fact that the 77^-orthogonal 
projection is stable in 77^(72), i.e., < |u|rri(n), 

11(1 - n;c)w||L2(0) 

K-1 

< \\{I -nK)u\\L2(n) + 

m—k 

< I e + ^ j |u|/ri(n) =: 

\ m—k / 


As hm = we obtain using the summation formula 

for the geometric series that T' < {e-\-2C]l‘^hk)\u\H':(Q,) 
for e —> 0 the desired result. ■ 

The statement © is well-known for all standard dis¬ 
cretizations. However, there might not be a good estimate for 
Ca- So, we are interested in the results by this lemma. The 
estimate ( [T6| ) can be treated using LFA. We can rewrite ( [T6| ) 
in matrix-vector notation as follows: 


W-PkK^'HKk+i)uk+i\\m^, < CAhi\\uk+^\^,, 


where Mk := := {J^ipk,j{x)ipk,i{x)dx)f^j^^ is 

the mass matrix. Here, the upper bound is obtained using the 
matrix norm: 


r^i/2 _ 


1 

hk 


- PkK-^PlKk+,)K-li^ 


Using the definition of the Euclidean norm and the fact that 
{I-PkK^^P^Kk+if = {I - PkK^^ P^ Kk+i), we obtain 


1 


Ca = ^p{ Mk+i{I - PkK 


-1 dT 


PfKk+,)K-l,). 


Qk + l- = 


Here, again, the spectral radius can be determined using the 
symbol 


Ca 


sup 

ee[o,2T7) 



where 


gYiih) ■-= (ww) ', 

cYiid) ■■= [i-K{e) [Kki0)y' Pkiori^iio)^ - 


As we have mentioned above, we are interested in comput¬ 
ing Ca for different discretizations. The details can be found 
in an accompanying Mathematica notebook, which is available 
in the wet]^ the main ideas will be given in the following three 
subsections. 


A. The Courant element 

The symbols JCk+i{d) and Pk{0) for the Courant element 
have already been determined in the last section. The mass 
matrix Mk has also a tridiagonal form. The symbol can be 
computed completely analogous as for the stiffness matrix: 


Mk+i{e) 


M^(0) _ 

+ tt) 


^ http://www.risc.jku.at/people/vpillwei/sLFA/ 


‘http://www.numa.uni-linz.ac.at/~stefant/J3362/slfa/ 










Fig. 6. Basis functions for the P^-spline discretization 


where 

^W = ^(e"'‘ + 4 + e«-). 
Based on this, we can derive 


Gk+i{0) — 


2 + COS0 —2 + COS0 \ 

—2 — cos 9 2 — cos 9 J 


The eigenvalues of Gk+i((^) 0 and 4. As this is already 

independent of 9, we immediately obtain that for the Courant 
element Ca = 5 is satisfied. 


B. A -spline discretization 

We can set up the same framework also for other discretiza¬ 
tions, like the discretization with splines. Here, assume that I 4 
is the space of all continuously differentiable functions, which 
are piecewise polynomials of degree 2. One possible basis for 
Vk is the basis of B-splines: 

forxi_i<x<Xi 
I - -t^{2x - Xi- Xi+if for Xi<x < Xi+i 
- Xi+2) for ^i+l <X< Xi+2 

0 otherwise, 

where Xi = ihk, see Fig. for a visualization of such a basis 
function. 

For the B-splines, we can again compute the integrals that 
are necessary to set up the mass matrix M^. As the support of 
the B-splines is larger than the support of the basis functions 
of the Courant element, we obtain a band matrix with a 
bandwidth of 5, with rrii^i = -^hk, 

Wi,i ±2 = Also for this case, we can determine the 

symbol 

^ (e-2‘® -f 26e-‘® -f 66 -f 26e‘® -f e^^®) . 

We can set up the the stiffness matrix and its symbol in 
a completely analogous way and obtain 

Kk{9) = ^ (-e-2'® - 2e-‘® + 6- 2e‘® - e^*® 

For setting up the symbol of the prolongation operator Pk-i, 
it is sufficient to solve again the equations ( [TT] i and ( [T2] i. For 
details, we refer to the Mathematica notebook. The overall 
symbol Qk+i{9) is again just obtained by multiplying the 
individual symbols. The eigenvalues of are 0 and 

—51-I-14cos(20)-I-cos(40) 

40(—2 -I- cos(0))(2 -I- cos(0))(2 -f cos(20)) ’ 


Fig. 7. Basis functions of the first kind of the -discretization 
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Fig. 8. Basis functions of the second kind of the P^-discretization 


This second eigenvalue can be rewritten using the replacement 
COS0 —>■ c as rational function, where the terms cos(20) 
and cos(40) are treated using the corresponding Chebyshev 
polynomials. Here we obtain - using CAD - that | is the 
largest value taken by ( fTS] ), so we obtain Ca = §• 

C. A standard P^-discretization 

Besides the spline functions, there is another possibility of 
setting up a discretization based on polynomials of degree 
2, which is even more popular in finite elements; we define 
14 to be the space of continuous functions that are piecewise 
polynomials of degree 2. Here, we can introduce a nodal basis, 
i.e., a basis where each basis function is associated to node 
(this basis function takes the value 1 on that node and the value 
0 on all other nodes). Here, the nodes are allocated on the ends 
of the intervals (as for the Courant element) and, additionally, 
on the midpoints of the elements. Here, we have two types of 
basis functions, cf. Fig. [7] and Fig. for visualizations. 

Because there are two types of elements, the mass matrix 
has alternating coefficients, see the Mathematica notebook for 


details; 
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For determining the symbol of Mk, we rewrite Mk as a sum of 
a band-matrix and of a residual matrix with alternating signs: 

l^k = , 

where Ak = is a band matrix with 

ai,i±i = j^hk and ai,i ±2 = --^hk and Bk is a matrix with 


















alternating coefficients: 
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V J 

Based on this decomposition, we can find the symbol. The 
symbol of Ak is obviously just 

Ak{e) = + 4e-i®‘ + 24 + 4e®' - e^®'). 

60 

The symbol corresponding to Bk is determined as follows: 

= (2 + e^^^+£^^,(0 + 7r ). 

Sl{e)-.= 

So, we obtain 

Mfe4(0) = Ife (0)4(0) + (0)4(0 + 4 

and, as 0 + 27r ~ 0 , also 

Mk^^{6 + tt) = Bk{9 + 7r)4(0) + Afe (0 + + tt). 

This shows, that Aik does not preserve a one dimensional 
linear span anymore, but a two-dimensional span, spanned by 
0^(0) and (j)^{d+TT). This is similar to the coarse-grid operator 
in the last section and in the last two subsections. So, the 
symbol is a representation of Mk with respect to the basis 
formed by these two vectors: 


KiO) = 


^k{0) _^( 0 ) \ 

Sfe(0 + 7r) AkiO + ir) j 


symbol Pk-i is a 2 x 4-matrix, for details we refer to the 
Mathematica notebook. Based on the symbols of the individual 
components, we can again compute Qk+i{9), the symbol of 
the overall operator. The eigenvalues of this matrix are 0,0, 
and so we obtain Ca = j|j- 

So, we have seen that the constant Ca takes the value 
I for the Courant element, the value | for the P^-spline 
discretization and ^ for the standard discretization. 

This indicates that the standard P^ discretization has the 
best approximation properties. However, the standard P^ dis¬ 
cretization needs two degrees of freedom per element, while 
the other two discretizations need, each, one degree of freedom 
per element. By dehning hk to be the distance between two 
nodes, i.e., hk = \hk for the standard -discretization and 
hk = hk for the other two discretizations, we can redehne the 
approximation error estimate as follows: 

||u — nfeu||4(n) < ^Ahl\uffji^Qy 

Here, we obtain Ca = | for the Courant element and Ca = | 
for both of the quadratic discretizations. 

As we have already mentioned, an extension to two dimen¬ 
sions is possible, however the terms get much more compli¬ 
cated. We refer to the complementary material, where we made 
an attempt to generalize the analysis to two dimensions. 

V. Concluding remarks 

We have seen that the terms that are constructed using LFA 
can be treated well using symbolic computation, particularly 
using CAD. Moreover, we have seen that the method of 
LFA can be applied in a wide range of problems. Besides 
is application to multigrid solvers, which is well studied in 
literature, cf. m, na, Eoi, LFA can be applied to other 
problems occurring in numerical analysis, like the computation 
of approximation error estimates. 
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